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Solar Flare Physics 

We have continued our previous efforts in studies of fourier imaging methods 
applied to hard X-ray flares. We have performed physical and theoretical 
analysis of rotating collimator grids submitted to GSFC for the High En- 
ergy Solar Spectroscopic Imager (HESSI). We have produced simulation 
algorithms which are currently being used to test imaging software and 
hardware for HESSI. We have developed Maximum-Entropy, Maximum- 
Likelihood, and “CLEAN” methods for reconstructing HESSI images from 
count-rate profiles. This work is expected to continue through the launch 
of HESSI in July, 2000. 

Section 1 of the following pages shows a poster presentation “Image Recon- 
struction from HESSI Photon Lists” at the Solar Physics Division Meeting, 
June 1998; Section 2 shows the text and viewgraphs prepared for “Imaging 
Simulations” at HESSI’s Preliminary Design Review on July 30, 1998 . 



SECTION I 

Image Reconstruction from HESSI Photon Lists 
Poster presentation at the 
June 1998 Solar Physics Division Meeting 
Boston, Massachusetts 



IMAGE RECONSTRUCTION FROM HESSI 


PHOTON LISTS 

1. PHOTON LISTS 



O WHY MAKE HESSI SCORES? (.ps) 

O 900-PHOTON SCORE (.ps) 

O PHOTON LIST RECIPE (.ps) 

O ALGORITHM DIAGRAM (.ps) 

O BINNED PHOTON PROFILE, COLLIMATORS 1-4 
O BINNED PHOTON PROFILE, COLLIMATORS 6-9 

2. BACK PROJECTIONS 

O Unbinned Back Projection (Caption) (.ps) 

O Unbinned Back Projection maps 0-8 (.ps) 

O Back Projection HOW-TO (.ps) 

O Back-Projection Artifacts (.ps) 

3. CHI SQUARE AND ALL THAT 

O Comparing the Model with the Score (.ps) 

O CHI-SQ for Sparse Photons (.ps) 

O Unbinned Back Projection maps 0-8 (.ps) 

4. MAXIMUM ENTROPY 

O Max Ent Algorithm (.ps) 

O HOW-TO MEM (.ps) 

O The Gradient Map (.ps) 

CONCLUSIONS 

This is mainly a progress report and has no final conclusions, but we are able to point out the directions 
for future work: 



• Photon lists must be developed to include: 

O More complex and realistic sources 
O Time variability 
O Internal shadowing and vignetting 
O Grid irregularities 
O Pulse height spectra 
O Aspect Phase Information 

• Maximum Entropy Programs must be made: 

O More robust 
O Faster 
O More efficient 

• Other Reconstruction Algorithms are needed: 

O CLEAN 

O Maximum Likelyhood 
O Signal Processing 


Last updated JUNE 2, 1998. 




PHOTON LIST SIMULATION 


A photon list is a table of arrival times, detector numbers, and 
pulse heights for time- tagging telescope such as HESSI. A simu- 
lated photon list (also called a “score”) can be used for: 

• Testing HESSI hardware during development stages 

o Simulated throughput with a “pulser stimulator” 
o Simulation of dead-time effects, pulse pileup, etc. 

• Testing HESSI software 

o Proving that mapping routines work 
o Identifying needs/uses/problems in utilities 
o Estimating times and memory required for analysis 

• Providing a realistic, simple, model of HESSI output 

o for scientific users 

o and the general public (educational-outreach) 

We have provided a simple score file in FITS format, along with 
an IDL program to read it and create a simple image. The file 
”scorelK0-4.fits” is available on the hesperia ftp site, and can be 
accessed using anonymous ftp: 
ftp hesperia.gsfc.nasa.gov 
cd pro 
binary 

get scorelK0-4.fits 
get score2bp.pro 
bye 

Alternatively, you can use a browser: URL = http://hesperia. 
gsfc.nasa.gov/hessidoc/pro/, and download the aforementioned 
.fits and .pro files. The FITS file is displayed with the IDL pro- 
gram score2bp.pro, which back- projects the data onto 64 x 64 
arrays. In IDL just type the command “.run score2bp”, and you 
will be prompted for the name of the score file. 

The following figure shows a 900-photon subset of the score file. 
Each “staff” shows 9 one-second timelines, with spatial frequency 
increasing up (collimators 1-9) and time going to the right. The 
photon arrival times from a point-source burst have been modu- 
lated by an analytical version of HESSI grids, and images can be 
constructed from the arrival times alone, as shown below in this 
paper. 









RECIPE FOR A PHOTON LIST 


1. Select a flare model in space and time. 

2. Create a pseudo-random time sequence tj within the desired 
time interval. 

3. Compute the rotation angle sequence 4>\ for one collimator. 

4. Calculate the HESSI modulation patterns for each (f >\ . 

5. Generate a sequence of relative probabilities Sj from the scalar 
product of the model and the modulation pattern at each <f) j. 

6. Create a uniform random number sequence Ri of the same 
length, scaled to the range [0, max(p;)]. 

7. Accept each time/angle value associated with those of Si which 
are smaller than the corresponding Rj. (I.e, ’’modulate” the pho- 
tons.) 

8. Append the accepted times, along with their collimator num- 
bers, to the photon list. 

9. Repeat until done. 



PHOTON LIST ALGORITHM 
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BINNED PHOTON LISTS 


COLLIMATOR 1: TBINS-3600, FILE^'SCOREI K_0-4.FITS" 
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COLLIMATOR 2: TBINS=21D0> FILE="SCORE1 K_0—4.FITS M 
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COLLIMATOR 3: TBINS=1215> F1LE="SC0RE1 K_Q-4.FITS" 
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COLLIMATOR 4: TBINS-700, F1LE= M SC0RE1 K_0-4.FITS" 
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BINNED PHOTON LISTS 


COLLIMATOR 6: TBINS=240, FILE="SCORE1 K_0-4,FITS M 
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COLLIMATOR 7: TBINS=135, F1LE="SC0RE1 K_0— 4»FITS" 
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UNBINNED BACK PROJECTIONS 

The following figures show full-disk (note Hmb circle) and small- 
field maps made directly from the score. Reading from left to 
right, these show projections from the coarsest 4 collimators and 
their sum (map 6). A box in the latter map shows the field 
mapped by the finest four collimators. These are shown indi- 
vidually for collimators 1-3, followed by their sum with the 4th 
collimator in the last map. 





















Artifacts in Back-Projection Maps 

• It has been shown that “Back -Projection” produces 
the statistically best estimate of the source position. 

• The Back-Projection map, of course, is not the “clean- 
est” or simplest possible map representing the true 
source. 

• If the collimator phase is constant during rotation, the 
Back-Projection map shows a ficticious source at the 
spin axis position due to the lack of modulation at that 
point. 

• The collimator phase will almost certainly not be con- 
stant for HESSI, so the spin axis will be washed out. 
Thus Back-Projection maps may be better in this re- 
spect than shown here. 

• There is a fictitious “mirror” source seen in Back- 
Projection maps at a position symmetric to the true 
source relative to the spin axis. The sign of the mirror 
source depends on the collimator phase, and its vari- 
ations will also wash it out with respect to the true 
source. 

• The rings around the source in Back-Projection maps 
have radii and spacing that depend on the collimator 
pitch. Adding the maps from several collimators there- 
fore tends to smooth out the rings, while reinforcing 
the true source. 

• The above being said, there is no substitute for a truly 
“clean” map, and Back-Projection maps will provide 
only the first look, after which reconstruction will take 


over. 



COMPARING MODELS WITH THE SCORE 


• The score can be thought of as a time profile with 
all zeroes and ones {n«}, the zeroes being as signif- 
icant as the ones. 

• In general, the model count rate profile {e«} will 
be non integral and non zero. 

• The x a statistic can be used as a measure of fit of 
{e«} to {n,} but it will not follow the x a distribution 
if n,- < 10. 

• Binning can mitigate this problem, but can result 
in loss of real temporal structure. 

• Bin sizes should be determined by the collimator 
geometry, the burst location and the rotation rate, 
not the number of counts. 

• Alternates exist to the x a statistic, such as the cor- 
relation coefficient or the C statistic. 

• Fbr low count rates we use the C = -2 log(P) statis- 
tic (Cash 1979, Ap.J. 228, 939), where: 


*=n 




• This is the probability that a Poisson-distributed 
sample {n,-} arises from the measured (or modeled) 
sample {e<}. 



The Chi-squared Statistic in the Poisson Limit 

In comparing a model e; of the observed photon count rate nj, it 
is common to employ the statistic: 

N t \ 2 

S _ ( n i - ej) 2 

. ej 

i=l 1 

In the limit of large values of n i? S follows the y 2 distribution, but 
it does not when the values of the countrate nj are all ones and 
zeroes, or if many values of ni are very small. Since we expect 
that many gamma-ray events observed by HESSI will have only 
a few photons, whose distribution follows Poisson statistics, we 
must consider the best way to find confidence limits of models. 
Following Cash (1979), we adopt the statistic: 

C = — 2(nj ln(ei/n 5 ) - e 5 + ni). 

It is easily shown that C -» S as n; -)• oo. In practice, C « S if 
ni > 10, and like S, C vanishes if ej = nj. 

For the case where there are p parameters and p-q constraints, 
one obtains confidence limits by letting all p parameters vary, 
finding the mimimum (C p ) min , and by varying p-q parameters to 
find the minimum (C p _ q ) min . Then the statistic 

AC = (C p _q) m ; n — (Cp) min 

obeys the y 2 probability distribution with q degrees of freedom. 

In particular, if one wants to find confidence limits for a particular 
pixel of a Maximum Entropy map, it is necessary to perform two 
minimizat ions : 

1. Compute the quantity (C p ) min from the {ej} found by the 
Maximum Entropy search through parameter (pixel) space. 

2. Select one of the pixels to test for goodness of fit, then de- 
termine (C p _!) min by searching for the minimum of C using 
MEM for all but that one pixel, whose flux is held constant. 

Finally compute AC, and determine the confidence level from the 
y 2 distribution with one degree of freedom. By Wilk’s theorem, 
this is valid as long as the total number of photons is large (>> 10), 
even if individual count rates are all zeroes and ones. 



MAXIMUM ENTROPY ALGORITHM 











MAXIMUM ENTROPY ALGORITHM 

1. Start with a blank map B, and a small positive value of pa- 
rameter A. 

2. From B, compute the model modulation profile {ej}. 

3. Calculate the C (or x 2 ) statistic measuring the difference be- 
tween the observed modulation profile {nj} and {ei}. 

4* Prom the modulation pattern, {n;} and {ei}, compute the gra- 
dient of C (or x 2 ) with respect to the free parameters (the map 
pixels). This is the “gradient map”. 

5. Using the fixed-point method (Sakao, 1992), find the map B 
which yields a maximum of the difference of the entropy and A 
times the gradient map. 

6. If C (or x 2 ) has not fallen below a pre-determined limit, in- 
crease A by a few percent and repeat from step 2. 

7. B is the “Maximum Entropy” map. 



GRADIENT MAPS 


One of the most important parts of the Maximum Entropy algorithm is the 
computation of the grad y 2 map. This is the gradient of the “goodness- 
of-fit” function y 2 or, in our case, C with respect to the map pixels. It is 
a reasonably good representation of the true map, and therefore deserves 
mention. 

The gradient of C with respect to the model modulation profile is very 
simple: 


dC_ 

dei 


1 - n?/e? 


Its derivitatives with respect to the map pixels /; are then given by the 
chain rule: 


dC _ ^ dC de k 
dfi ““ dek dfi 

The last derivatives on the right hand side are simply the modulation pat- 
terns, because the e* are simply the inner product of the modulation pat- 
terns with the model map /;. 

The figure below shows the first gradient map of the Maximum Entropy 
for the same case as the back-projections above. Note the absence of the 
spurious spin-axis source. 



SECTION II 
Imaging Simulations 
Presentation at 

HESSI’s Preliminary Design Review 
July 30, 1998 
College Park, Maryland 



IMAGING SIMULATIONS 


WHAT WE HAVE DONE 


Modulation Patterns 

• Analogous to HXT, but for HESSI will vary according to application. 

O HXT has only 64 modulation patterns, one per subcollimator. 

O HESSI will have hundreds to thousands. 

O FOV size and location are arbitrary in HESSI. 

O Center of rotation not constant, so pattern phases will vary with time. 

O The number of patterns equals the number of time bins. 

• Essential to simulations and imaging software 

O Initial maps are weighted sums of the modulation patterns. 

O Even pure Fourier methods have to return to image space and add patterns together. 
O Modulation patterns exploit HXT, HEIDI heritage. 

• Patterns derived from empirical measurements 

O Presently are based on theoretical grid models. 

O By April 1999 will be derived from OGCF and XGCF measurements. 

O Grid imperfections will be incorporated using phase and amplitude matrices. 

• In general, patterns need to be constructed quickly. 

O They form an "inner loop" in most reconstruction algorithms. 

O Impractical to rotate fixed patterns— slower than on-the-fly computation. 

O Every pattern must be phase shifted using the aspect solution. 

O Memory/disk constraints may mandate "on-the-fly" computation. 

■ 1000 (say) 128 x 128 patterns require 65 MB of memory; 

■ A slow computation usually beats all but fastest disks; 

■ On-the-fly computation provides more freedom of action. 

• Recent breakthrough: Order of magnitude speedup of computation: 

O Method involves copying from a 1 -d profile of the subcollimator profile 
O A graphic explanation 

O Makes it possible to back-project scores in < 4 s per collimator. 

• Free parameters of the modulation pattern 

O Slit/Slat ratio, a function of distance from spin axis, 

■ Internal shadowing important — can be a 10-20% effect 

O Harmonic content: Fundamental and any subset of harmonics 1 to ao, 

O Arbitrary phase offsets given by aspect solution, 

O Phase and amplitude correction matrices. 


» 




This figure shows the 64 modulation patterns of the Hard X-ray Telescope on the Yohkoh satellite. Each 
square of black and white bars is the visibility pattern (Fourier pattern) for a single sub-collimator. The 
finest patterns have a spatial resolution of 8". 


This figure shows a small subset of the modulation patterns for the HESSI Telescope. Each column of 
squares is a subset of the visibility pattern (Fourier pattern) for a single sub-collimator during one 
telescope rotation. The counts are binned for 0.5-40 millisec to provide sufficient statistics for one 
Fourier component. The finest patterns, which are binned for about 2.5 ms, have a spatial resolution of 
2.3". Only the first 12 of 1620 are shown. The coarser collimators are binned for progressively longer 


intervals in proportion to their spatial resolution (2.3, 4.0, 6.9, 12.0, 20.7, 35.9, 62.1, 108, and 186 arc 
sec, resectively). The number of modulation patterns depend on the degree of binning chosen, but in this 
case, if all were shown, the number of matems in each column would be 1620, 935, 540, 312, 180, 104, 
60, 35, and 20, respectively. 



LOOK-UP TABLB ARRAYx 


MODULATION 
PATTERN 
ARRAY 


Bach pixel of the modulation pattern array is created by reading a 
value out of the lookup table F(s), which is obtained from the 
collimator pre-flight measurements. F(s) is shown above as the 
curve to the left of the thick line OS, -which lies at an angle 9 
measured CCW from the x axis 

General points in the array are taken by projecting 

the perpendicular to the line OS and reading F(S). This method 

replaces a cosine by a lookup, which is an order-of-magnitude 



Shadowed *>lit/'pit.:h rati:* 


Effect of shadowing on s lit / pit-; h ratios 
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Infinite photon case ( > 10000 events/coll/s ) 

• Maximum Entropy 

O Poisson noise was added to the model maps before reconstruction. 

O Gives nice results BUT... 

O Has NOT been tested using a score with realistic statistics. 

O X 2 statistic was not used in the earliest mapping (pre 1997) 

O Later (1997) X 2 statistic was used, (example) 

• Max Likelihood (Richardson-Lucy) —X 2 statistic was not used 

O To get convergence, the number of data points and map pixels must be approximately equal. 
O Not appropriate for low countrate regime. 

O Nevertheless, non-binning algorithm shows promise. 

• CLEAN 

O X 2 statistic was not used 

O But it can be incorporated to provide termination test. 

• Analogy with radio astronomy methods is good, BUT has its limits: 

O Radio techniques are not designed for Poisson-distributed photons. 

O In those routines which use statistical tests, a 2 ± count 

• All working programs will be revised to incorporate event lists so that they can work with real 
HESSI data. 
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MAXIMUM LIKELIHOOD RECONSTRUCTION 


Given a model omul; mix: time profile {Pi}, raid an oliwrwid count-rate profile {JVi}, wo 
wiili to determine the likelihood of tl s: otmerved jmjfile. (The N+ an: all integon anHumed 
to he Paumciudiittiihtitod, >izid the A an: not nootMuariiy intigeni, hut they an: all pmiifive.) 
From PoiHiKJu Hhathtira, the probability of getting N aouutH when t: an: expected in: 


p{N | D) = 


vxp(-D)D N 

JV! 


raid mo the joint likelihood L of getting the ohiicrvod countw N(i) iu the time hin i, given 
the: expected couutit D(i) i* given tjy: 


In L = ]T JY(i) in D(i) - D(i) - lnN[i)\ 
i 

The: maximum lilociihtxxl map wolutLon ocean when all jjarhial derivative** with nxpcct to 
the map pixels {Oi} winiak: 

din L _ fltn L d&(j) _ 

dO(f}-^ 6 D[j)dOU)- V 


The: quantity P(i j) in the modulatidn pattern, where the index j in the pixel number and 
the: Hum in over the time indoc i. The Rixhradxoii-Lucy iteration fix' tlx: Holutiitm xm: 


0^U) = 0(i)J2 r ^3)W: i 

t ^ ' i 


I 
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Finite photon case ( < 1000 events/coll/s ) 

• Scores 

O Simulations of the raw data structure— time, subcollimator, pulseheight 
O Important for testing software 

O And for running a "pulser stimulator" through the flight electronics. 

O Very simple Monte-Carlo program 
O Photon list algorithm diagram 

O See hessidoc ftp site: hesperia.gsfc.nasa.gov/hessidoc/pro. 

O Or WWW site: http://hesperia.gsfc.nasa.gov/~schmahl 
O Binned Photon Lists 

• Back projection— the basic imager 

O Back projection: How it works 
O Quick-look maps of flares 

O Useful to localize flare center for subsequent reconstruction 
O Can be used as "dirty maps" for "cleaning" 

• Reconstruction 

O Maximum Entropy, problems 
O Maximum Likelyhood, possibilities 

O In all methods, must define a count-rate based goodness of fit 

• Confidence limits: 

O Must define the appropriate statistic for count-rate profiles and maps. 

O X 2 Statistic not X 2 distributed 
O Webster Cash (Ap.J. 1979, 228, 939) showed the way. 

O Likelihood Function 

■ L=Tte c C N /N! 

■ Recipe for confidence testing of a pixel: 

1 . Run optimization/reconstruction using all map pixels as free parameters => L 

L ° 

2. Run optimization/reconstruction holding a map pixels at some test level => L 
L l 

3. Then In (Lj / L Q ) is X 2 -distributed. 

■ Details in SPD poster page 


» 
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WHAT WE PLAN TO DO 

Who is doing what? 

• More complex scores (Schmahl) 

O Exponentially increasing source 
O Multiple, extended sources 
O Incorporate simulated grid errors 

• MaxEnt, Max Likelihood, CLEAN (Schmahl) 

• Forward Model (Hurford) 

• PIXONS (McTieman?) 

• Fourier (Vilmer?) 

• Self-Calibration (Aschwanden) 

• Full-sun mapping (Kosugi?) 

O Needed for lowest energies 

O Big problems with over resolution for extended features 
O Internal shadowing varies across disk— How to incorporste it? 

• Spatial-Spectral 

O HXT sometimes shows dips in bremsstrahlung spectra (impossible) 
O Points out the need for spectral/spatial reconstruction 
O Can use penalty function (ala’ Gary 1996) or 3-D entropy 


E.J.S. July, 1998 


